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Abstract. We explore the regime of "superfast" reactivity that has been predicted 
to occur in turbulent flow in the presence of potential disorder. Computer simulation 
studies confirm qualitative features of the previous renormalization group predictions, 
which were based on a static model of turbulence. New renormalization group 
calculations for a more realistic, dynamic model of turbulence show that the superfast 
regime persists. This regime, with concentration decay exponents greater than that for 
a well-mixed reaction, appears to be a general result of the interplay among non-linear 
reaction kinetics, turbulent transport, and local trapping by potential disorder. 
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1. Introduction 

Reactive turbulent flow is important in a variety of natural processes, ranging 
from the production of smog in the atmosphere jlj to the feeding habits of certain 
oceanic creatures^. The mixing effects of turbulence are put to use in a variety of 
engineering applications, including combustion reactors, fluid catalytic cracking units, 
and polymerization reactors. The behavior of reactive turbulent flow is usually analyzed, 
numerically or analytically, with the continuum reaction/transport equations [I], |3[]. 
These mean-field equations fail, however, in two dimensions at low densities of reactants. 
This is understood in a general way, since the upper critical dimension for bimolecular 
reactions, even in the absence of turbulence, is two J|, |5], ||, [7j. 

A renormalization group treatment of a model of reactive turbulent flow has recently 
predicted a regime of "superfast" reactivity for the v4+v4 — > reaction in two dimensions 
||. The superfast reactivity occurs when a certain degree of quenched potential disorder 
is added to the turbulent system. In the superfast regime, the reactant concentration is 
predicted to decay more quickly than in a well-mixed reaction: i.e. c(t) ~ a/t 1+x at long 
times, with x > 0. This rapid decay occurs due to a subtle combination of the effects of 
turbulent transport, trapping by disorder, and reaction. It shows up in the analytical 
calculations as a renormalization of the effective reaction rate. It was anticipated that 
this result may have technological implications for certain thin-film reactors or fluidized 
MEMS devices. 

The renormalization group calculations were carried out with a statistical model of 
turbulence, in the same spirit as in treatments of the turbulent transport of passive 
scalars g, H, [TTJ . A review of this approach to the dynamics of a passive scalar 



can be found in [pL2| . The statistics of the turbulent velocity field in the model were 
chosen to reproduce the correct transport properties at long times. There was freedom 
to choose either a time-dependent or time-independent velocity field, and the time- 
independent choice was made. This choice corresponds to G. I. Taylor's hypothesis of 
frozen turbulence. 

In this paper, we perform further studies to confirm the existence of the superfast 
reaction regime. In Section 2, we describe a method for simulating reactive turbulent 
flow with random, time-independent velocity fields. In Section 3, we present and 
discuss our computational results. The existence of the superfast regime is verified, 
and additional features arising from higher-loop contributions in the renormalization 
group are found. In Section 4, we present a renormalization group treatment of reactive 
turbulent flow with a more realistic, time- dependent model of the velocity field. When 
the appropriate type of potential disorder is included, the superfast regime is found to 
persist. We conclude in Section 5. 



'Superfast" Reaction in Turbulent Flow with Potential Disorder 



3 



2. Simulation Methodology 

We consider the reaction 

A + A^V (1) 

occurring in a model of a two-dimensional, turbulent fluid. In the absence of reaction, the 
A particles are advected by the fluid streamlines. This chaotic motion is superimposed 
upon the natural Brownian motion of the particles. In addition, the particles experience 
a force due to the spatially-varying, random potential. In the presence of reaction, two 
A particles react with the conventional reaction rate constant A when in contact. 

We consider this reaction to proceed on a square N x N lattice, both for 
computational convenience and to allow a direct comparison with the field-theoretic 
results. How the configuration of reactants on the lattice changes with time is specified 
by a master equation. The master equation defines transition rates for all of the possible 
changes to a given configuration of reactants on the lattice: 

dP{{ ^ ht) = J^fo + l)P(...,n i - l,n 3 + 1, ...,*) - r^P] 

ij 

+^ E [fa + 2 )fa + 1 ) P (-' ^ + 2 ' "■'*) - n *fa ~ X ) p ] ' ( 2 ) 

i 

where n, is the number of A particles on lattice site i, and h is the lattice spacing. The 
rate of hopping from lattice site i to lattice site j is given by t^ 1 , which will be defined 
below. The summation over % is over all sites on the lattice, and the summation over j 
is over all nearest neighbors of site i. The particles are initially placed at random on the 
lattice, with average density uq. The initial concentration at any given site is, therefore, 
a Poisson random number: 

Pin,) = I^e-'^ 2 

(m/h 2 ) = n . (3) 

The turbulence and potential disorder enter this master equation through the 
hopping rates t~ 3 - . These rates are chosen so that the master equation reduces to 
the conventional transport equation when there is no reaction: 

— = DV 2 c + (3DV ■ [cVu - cV x 0] , (4) 

where c(r,t) is the concentration of the A particles, D is the diffusion constant, 
(3 = 1/(/cbT) is the inverse temperature, u(r) is the random potential, and <p(r) is 
proportional to the random stream function. A particularly simple form for the transfer 
rate from position % — r to the position j = r + Ar is 

^^[l-^.(V W -Vx^)/2] . (5) 

When equation (^) is used in equation (Q), one finds via a Taylor series expansion in h 
that the conventional transport equation (01) is reproduced. This form of the transfer 
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rates (H) is the one conventionally used in theoretical calculations || . In our simulations 
we choose, instead, to use the form 

= ^ exp{{3( Ul - Uj )/2} exp{/3(</y - 0,)/2} . (6) 

Here j and j' are nearest neighbors of site i. The identity of site j' is derived from site 
j by a counter-clockwise rotation of tt/2 about site i. This form leads to a transition 
rate in the simulation that is bounded and non-negative, in contrast to the simpler form 
of equation (||). A Taylor series expansion of equation @ leads directly to the simpler 
equation (|5|). Since both of these forms are equivalent in the limit of a small lattice 
spacing, they are expected to lead to identical scaling in the long-time regime. 

Following the usual statistical approach to turbulence, we take the stream function 
and potential to be Gaussian random variables with correlation functions given by 

Xuu{k) = ^ , (7) 

where the Fourier transform of a function in two dimensions is given by f(k) = 
j d 2 rf(r) exp(ifc • r). The random stream function and potential are generated first 
in Fourier space, where the values at different wave vectors are independent Gaussian 



random variables, and then converted to real space with a fast Fourier transform fT3 
Periodic boundary conditions are used on the lattice. The quality of the random number 
generator is of importance when simulating systems with long-ranged correlations fT3|j . 
We use a sum of three linear congruential generators method |T4J . 

Isotropic, fully-developed turbulence is modeled by y = 8/3 and 7 = 0. Developing 
turbulence is modeled by smaller, positive values of y. Potential disorder of the form we 
consider here could arise in a reaction between ionic species confined to a thin film of 
fluid between spatially-addressable electrodes or media with quenched, charged disorder. 
The ionic species would move in the random potential generated by the electrodes or 
quenched disorder. The electrodes or quenched disorder could be devised so as to 
reproduce the correlation function Xuuij)- 

The master equation (Q) can be exactly solved by a Markov Poisson process ||15|| . 
In such a process, we consider the motion of discrete A particles on the lattice. For any 
given configuration of reactants at time t, there are 4 rij possible hopping events. 
Each occurs with a rate given by equation (|^). In addition, there are Y2i n i( n i ~ -0/2 
possible reaction events on the lattice. Each of these occurs with the rate r"* = \/h 2 . 

The Markov process is started by initially placing the reactants on the lattice 
at random with average density n . Each step of the Markov process consists of 
randomly picking one of the possible reaction or diffusion events and incrementing time 
appropriately. The probability of event a occurring, out of all the possible diffusion and 
reaction moves, is 
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After performing the chosen event, time is incremented by 

where £ is a random number uniformly distributed between and 1. The Markov chain 
is continued until zero or one reactants remain on the lattice, at which point no more 
reaction events will occur. 

This Markov process, when averaged over initial conditions and trials, exactly 
reproduces the predictions of the master equation. Averages over the statistics of the 
turbulence and potential disorder are taken by directly averaging over sufficiently many 
instances of the random fields. Since the asymptotic scaling regime is reached only for 
long times, we make two approximations to facilitate the computations. First, we set the 
reaction rate to infinity. That is, when a particle moves to a site with another particle 
already present, the reaction occurs immediately. Physically, we do not expect this to 
modify the long-time decay law, since at long times the reaction will be in a transport- 
limited regime. Indeed, the renormalization group prediction for the decay exponent 
is independent of the reaction rate || . The predicted prefactor also has a well-defined 
value in the infinite reaction rate limit. Second, we assume that each particle on the 
lattice is equally likely to undergo a hopping event: 

Pfmoving particle a) = — , (10) 

n 

where n = Xli n « i s the total number of particles on the lattice. After the particle to 
move is chosen, the probabilities of each of the four possible hopping events are derived 
from equation @. In particular, another uniform random number is generated for 
comparison against the four different hop probabilities: 

P(hopz) = -^ r . (11) 

After the chosen hop is performed, time is incremented by At = N 2 /(n . t^ 1 ). The 
uniform choice of particles to move and the approximate time incrementation are not 
expected to influence the long-time exponents. In particular, these approximations are 
exact in the long-time limit if the hopping rates of equation @ are used. 



3. Simulation Results and Discussion 



The renormalization group calculations make a prediction for the concentration decay 
exponent, the a in c(t) ~ at~ a [|J. The predictions depend on the properties of the 
turbulence and potential disorder through the parameters y, a, and 7. The predictions 
for several values of y are shown in figure [l]. The decay exponent initially rises with 
increasing potential disorder and eventually decreases. The maximum reaction rate of 
a = l + y/(6 — y) occurs for a = 37. The decay rate depends on a and 7 only through 
the combination 7/a. These predictions come from a one- loop renormalization group 
calculation, and they are strictly valid only for small y. In addition, the predictions 
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should be more accurate for small 7, because it is known that higher-loop corrections 
lead to a modification of the flow diagram for large 7. 

We find that our Poisson process efficiently solves the master equation for values 
of y near unity. For smaller values of y, the renormalization of the effective reaction 
rate occurs slowly. This means that the predicted asymptotic scaling occurs only at 
long times, times longer than we can reach in our simulation. For larger values of y, 
significant lattice effects occur, due to our choice of transition rates, equation @, and 
correlation functions, equation (0). 

Shown in figure ||] are the decay exponents observed in our simulation for y = 3/4. 
The decay exponents observed for y — 1 and y = 5/4 are shown in figures ||] and ||, 
respectively. Each of the data points in the figure is an average over three different 
runs on three different instances of the turbulence and potential disorder. The standard 
deviations estimated from these three runs are shown as the error bars. 

We made a few canonical choices for constants in our simulation. So as to reach 
the asymptotic scaling regime, we used 4096 x 4096 square lattices. Simulations on 
2048 x 2048 lattices gave similar results. Finite size effects should be most noticeable 
for small values of y, a, or 7, since in that regime the renormalization of the effective 
reaction rate is most slow. Indeed, finite size effects only appear to be present in figure 
|2] for small 7/0", where the slow renormalization of A(/) leads to an observed decay 
exponent that is below the correct asymptotic value of unity. While additional finite 
size effects could be present, the agreement at small 7/cr among curves for different 
cr but same 7/cr and y, and the agreement between the results for 2048 x 2048 and 
4096 x 4096 lattices, argues against this. We used a lattice spacing of h — 1, which 
we can enforce by a spatial rescaling, and a diffusion constant of D = 1, which we 
can enforce by a temporal rescaling. The long-time decay exponent is not affected by 
these rescalings. We collected the time-dependent concentration data on the lattice by 
binning the results into a histogram with a temporal bin width of St = 1. We continued 
each simulation until the concentration was so low, c ~ (const) /iV 2 , that lattice effects 
were obvious. We show in figure |5] the simulation results for a typical run. The average 
of the slopes determined in three such runs gives one of the data points in figures |2|44|. 

Since the renormalization group predictions depend only on 7/cr, we fixed a and 
varied 7 for each simulation at a fixed value of y. Under these conditions, the effective 
strength of turbulence is expected to flow to its fixed point value, cr — > a* « 2iry(27i/h) y 
||. In each of the figures ||-[| we show results for three different "bare" values of a: 
a = a\ a = a*/2, and cr = 3cr*/2. 

We use a lnc(i) versus hat plot to obtain the decay exponent. This plot has many 
more data at large hit, which leads to an over-weighting of the long-time regime. To 
counteract this effect, we use data exponentially spaced in time. We exclude both 
short-time data, which are not in the asymptotic regime, and long-time data, which 
show finite-size effects. Use of exponentially spaced data allows one to obtain a reliable 
estimate of the error in the exponent determined for a given realization of the disorder. 
We find that this statistical error is smaller than the systematic deviation that occurs 
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for each different realization of the disordered streamlines and potential. The error bars 
shown in the figures encompass both the statistical and the systematic errors. 

We see that there is a general agreement between the simulation results and the 
renormalization group predictions. In particular, the superfast reaction regime, a > 1, is 
observed. This was a dramatic analytical prediction, and confirmation of this regime is 
significant. The exponent is observed to reach roughly the maximum value of 1+y / (6—y) 
predicted by theory. 

The agreement between the simulation results and the one-loop renormalization 
group predictions is better for smaller y. For small y, the location of the peak in 
reactivity is trending towards the predicted value of 7/cr = 1/3. For finite values of y, 
the location of the peak is shifted to the left. Indeed, the whole curve is shifted to the 
left for finite y. One consequence of this is that the observed decay exponent decreases 
below unity for 7/cr close to, but less, than one, in contrast to the one-loop predictions. 

The agreement between the simulation results and the one-loop predictions is also 
better for small 7. For small 7/cr, the simulation results for different values of a fall 
on the same curve. This is the universal curve predicted at one-loop, albeit compressed 
to the left. For large 7, the simulation results depend on both 7 and a, not simply on 
the ratio 7/cr. This dependence is presumably due to corrections that would enter in a 
higher-loop calculation. 



4. Renormalization Group Treatment of Dynamic Turbulence 

The simulations, and the previous renormalization group calculations, were performed 
for a model of turbulence with streamlines random in space, but constant in time. For 
this model, both the analytical and computational studies predict a regime of superfast 
reaction in the presence of potential disorder. This regime was understood to occur 
as a general consequence of the interplay among turbulent transport, trapping by the 
random potential, and a non-linear dependence of the reaction rate on local reactant 
concentration. This regime was not believed to occur, for example, due to some subtle, 
artificial interaction between the quenched potential and quenched streamlines. 

Even so, it would be interesting to predict the existence of the superfast reaction 
regime for a dynamic model of turbulence. As we have seen, reliable simulations require 
large lattices and long times, even for quenched turbulence. Moreover, simulations were 
possible only for an intermediate range of y, due to the competing considerations of 
minimizing lattice effects and accessing the asymptotic regime. Renormalization group 
theory is, however, a viable approach to studying reactive turbulent flow in a dynamic 
model of turbulence. 



We, again, use a statistical model of turbulence JTTJ] . We now assume that the 
stream function is random in both space and time: 

x**{k, h - h) = - t 2 \- p Q{h - t 2 ) , (12) 

where G(t) is the Heavyside step function. A family of models for turbulence is generated 
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by varying p and y a . Isotropic turbulence is modeled when 3y a — 1p 
small p. We include a random potential, as before: 



Xuu(k) 



7 



k 2 +y~t 



8, at least for 



(13) 



We will find that the most interesting regime occurs for y 7 < y a , due to the weakening 
of the random stream function through decorrelation over time. 

As in previous studies ||, we map the master equation (0) onto a field theory 
using the coherent state representation 0, [L6[ . For this operation, we use the transition 



rates given by equation ([5]). The random and stream function are incorporated with 
the replica trick | 10[| , using N replicas of the original problem. The concentration of 
reactants at time t, averaged over the initial conditions, c(r,t), is given by 



c(r,t) 



lim (a(r, t)) 

7V->0 



(14) 



where the average is taken with respect to exp(— S), with 

S = V / d 2 r / dta a (r, t) [d t - DV 2 + <*(*)] a a {r, t) 

2a a {r, t)a 2 a {r, t) + a 2 (r, t)a 2 a (r, t) 



a=l 

N 

a=l 

(3 2 D 2 



dt 



d ra a (r, 0) 



j dt idt 2 j (27r) 2 5(fc 1 + k 2 + fc 3 + fe 4 ) 



xa ai (fei, tx)a ai (fc 2 , h)a a2 (fe 3 , t 2 )d a2 (fc 4 , t 2 ) 
x ki ■ (fei + k 2 )k 3 ■ (fe 3 + ki)x uu (\ki + k 2 \) 

+ki x k 2 k 3 x kiX H (\ki + fc 2 |,t 2 - *i) • (15) 

The notation J fc stands for j d 2 fc/(27r) 2 . The upper time limit in the action is required 
only to satisfy U > t. 

Using renormalization group theory, we can derive the flow equations for this 
model. We find that the turbulence contributes directly to the dynamical exponent, 
and therefore to the transport properties, but not directly to the effective reaction 
rate. The potential disorder contributes directly to both the dynamical exponent and 
the effective reaction rate. Finally, the dynamical exponent contributes indirectly to 
the effective reaction rate as well. We define g a = j3 2 aA 2p ~ y<T D P T(1 — p)/(47r) and 
g 7 = (3 2 jA~ y ~< /(47r), where A = 2n/h is the cutoff in Fourier space, and Y{x) is the 
standard Gamma function. The precise definition of g a depends on the regularization 
°f Xh at t = 0. At one loop, we find the flow equations to be 
dlnn = 
dl 
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3g 7 - 9a 



dl AttD 
d\n g 7 



dl 
d\n g a 



y„- pz- 2g a . (16) 



dl 

The dynamical exponent is given by 

z = 2 + g^-g a . (17) 

To achieve the superfast regime, we must set y 1 = y a — pz. Without this choice, 
the potential disorder will either dominate or be irrelevant. It is natural that the 
appropriate y 1 should be weaker than y a , since the effective turbulence strength is 
weakened by decorrelation effects over time. With this choice, we find that the ratio 
9j{l) / 9<r(l) remains constant under the renormalization group flows. The parameter 
g a (l), however, flows to the fixed point g* = y 7 /2. Moreover, we find that when 
9-yl 9a > 1/3, the effective reaction rate flows to a finite fixed point A* = 4:7iD(3g* — g*). 
When g 1 /g a < 1/3, the effective reaction rate renormalizes to zero. 

With the choice of y 7 = y a — pz, the flow equations are the same as for the static 
model of turbulence ||. The matching is the same, as well. In terms of these new 
variables, therefore, the predicted concentration dependence is the same. In particular, 
the superfast reaction regime is predicted to exist for this dynamic model of turbulence. 
The explicit, long-time concentration dependence for weak disorder is 

-2gl/(2+gl-g*) 



C(t) 



1 1 

+ 



[A7rD(g*-3g;) A 



1 / t 



t \t 



o 



(3s 7 < g a ) , (18) 



where Ao is the bare value of the reaction rate. The matching time, to, is given roughly 
by to ~ h 2 /(2D). For strong disorder, the concentration decays as 

<<*) ~ „_„,o_ ^A t) , (3^ 7 >^). (19) 



4irD(3g*-g*)t \t 



o 



At the location of the maximum reaction rate, there is a logarithmic correction to 
power-law decay: 

c ( t ) ~ Ht/to) t -y 7 /(6-y 7 ) f3 \ (2Q) 



5. Conclusion 



A regime of reaction rates faster than that for a well-mixed reaction was observed 
in simulations of a static model of turbulence plus potential disorder. Qualitative 
predictions of previous renormalization group studies of this model were reproduced. 
Interesting departures from the one-loop predictions, such as a compression of the decay 
exponent curve to the left, were observed. These features are expected to be reproduced 
in higher-loop analytical calculations. 
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Renormalization group calculations for a more realistic, dynamic model of 
turbulence show that the superfast regime persists. Indeed, these calculations suggest 
that our general understanding of the phenomenon of superfast reactivity is correct. The 
random potential tends to attract reactants to localized regions in space, and by doing 
so increases the reaction rate in a non-linear manner in these regions. The reactants, 
therefore, annihilate at a faster than average rate in these regions. The turbulent velocity 
fields transport reactants to these regions rapidly, although the transport is slowed by 
the ruggedness of the random potential landscape. The net effect of the trapping and 
the turbulent transport is superfast reactivity whenever the mean-square displacement 
is super-linear. 

The existence of a superfast reaction regime is of great interest from a practical point 
of view. Often, it is assumed in reactor design that the highest reaction rate is obtained 
for a well-mixed system [|17]. Our results show that under certain conditions, the highest 
reaction rate is achieved, instead, for an inhomogeneous system. This behavior has 
been observed experimentally. Indeed, an enhancement due to inhomogeneous reactant 
concentrations of the effective reaction rate between ions in a chaotically-mixed, two- 



dimensional, fluid with attractors has recently been observed [0 . The superfast regime 
may be relevant to a variety of reactors, including thin-film reactors, certain combustion 
reactors, and certain types of micro-electro-mechanical (MEMS) devices. 



Acknowledgments 

This research was supported by the National Science Foundation through grants CHE- 
9705165 and CTS-9702403. 

[1] R. W. Bilgcr, L. R. Saetran, and L. V. Krishnamoorthy, J. Fluid Mech. 233, 211 (1991). 

[2] T. Ki0rboe and E. Saiz, Mar. Ecol. Prog. Ser. 122, 135 (1995). 

[3] J. C. Hill, Ann. Rev. Fluid Mech. 8, 135 (1976). 

[4] L. Peliti, J. Phys. A: Math. Gen. 19, L365 (1986). 

[5] B. P. Lee and J. Cardy, J. Stat. Phys. 80, 971 (1995); 87, 951 (1997). 

[6] J.-M. Park and M. W. Deem, Phys. Rev. E 57, 3618 (1998). 

[7] M. W. Deem and J.-M. Park, Phys. Rev. E 57, 2681 (1998). 

[8] M. W. Deem and J.-M. Park, Phys. Rev. E 58, 3223 (1998). 

[9] D. S. Fisher, M. P. A. Fisher, and D. A. Huse, Phys. Rev. B 43, 130 (1991). 
[10] V. E. Kravtsov, I. V. Lerner, and V. I. Yudson, J. Phys. A: Math. Gen. 18, L703 (1985). 
[11] M. Avellaneda and A. J. Majda, Phys. Rev. Lett. 68, 3028 (1992). 
[12] J. P. Bouchaud and A. Georges, Phys. Rep. 195, 127 (1990). 
[13] V. Pham and M. W. Deem, J. Phys. A: Math. Gen. 31, 7235 (1998). 
[14] B. Wichmann and D. Hill, Byte 12, 127 (1987). 
[15] T. Elston and C. Doering, J. Stat. Phys. 83, 359 (1996). 
[16] B. P. Lee, J. Phys. A: Math. Gen. 27, 2633 (1994). 
[17] K. Jensen, Nature 393, 735 (1998). 

[18] O. Paireau and P. Tabeling, Phys. Rev. E 56, 2287 (1997). 



"Superfast" Reaction in Turbulent Flow with Potential Disorder 
Figure Captions 

1.30 



1.20 



11 



1.10 



1.00 



0.90 



0.80 




0.0 



0.5 



1.0 

y/a 



1.5 



2.0 



Figure 1. The decay exponent for the A + A — > reaction: c(t) ~ (const)t a . The 
correspond to y = 1 (solid), y = 5/4 (short-dashed), and y = 3/4 (long-dashed). 
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Figure 2. The decay exponent observed in the simulation for y = 3/4. The three 
curves correspond to a = a* (solid), a = 3<r* /2 (short-dashed), and a = cr*/2 (long 
dashed) . 




Figure 3. The decay exponent observed in the simulation for y = 1. The three curves 
correspond to a — a* (solid), a = 3a* / 2 (short-dashed), and a = a* /2 (long-dashed). 




Figure 4. The decay exponent observed in the simulation for y = 5/4. The three 
curves correspond to a = a* (solid), a = 3<r* /2 (short-dashed), and a = cr*/2 (long- 
dashed) . 




Figure 5. The concentration as a function of time (solid line) for the case y = 1, 
a = a*, and 7/cr = 0.2. Also shown is the fit to the data in the scaling regime (dashed 
line) . 



